Samenvatting

Om te laten zien dat GIS mogelijk is in R, gebruiken we wat open geodata van Den Haag, en maken een map.

Inleiding

Den Haag heeft een aantal kaarten openbaar gemaakt. Wij gebruiken de kaart van blokverwarming, energielabels (p5), en combineren deze met open satellietdata van NASA, en een open street map.

Data

Eerst downloaden en upzippen we de shapefiles van geoportaal-ddh.

# file download & unzip
download.file("http://geoportaal-ddh.opendata.arcgis.com/datasets/9cca53a6a4094a0b80e964b181a484ad_3.zip","blokverwarming.zip", mode="wb")
unzip("blokverwarming.zip")
download.file("http://geoportaal-ddh.opendata.arcgis.com/datasets/ed833f69fb82463ea63d37be00c356c3_0.zip", "energielabels.zip", mode="wb")
unzip("energielabels.zip")

# inputfileblokverwarming
f_blokverwarming <- "blokverwarming_DH.shp"
# inputfile enegielabels
f_energielabels <- "Energielabels_Postcode_5_niveau_Den_Haag_2016.shp"

Kaart Energielabels

We lezen de file in met readOGR van package rgdal. En daarna plotten we de polygons met wat kleurtjes, met behulp van base plot.

require(rgdal)
# de layers uit de shapefile
fl_energielabels <- ogrListLayers(f_energielabels) ## het is er slechts 1
# lees de layer in met readOGR
energielabels <- readOGR(f_energielabels,layer=fl_energielabels[1])
## OGR data source with driver: ESRI Shapefile 
## Source: "Energielabels_Postcode_5_niveau_Den_Haag_2016.shp", layer: "Energielabels_Postcode_5_niveau_Den_Haag_2016"
## with 844 features
## It has 10 fields
## Integer64 fields read as strings:  OBJECTID FREQUENCY labelCodeR MaxSimpTol
# plot de kaart
plot(energielabels,col=colorRampPalette(c("brown", "grey"))(5))

We willen natuurlijk graag laten zien wat de gemiddelde energielabels zijn per p5 gebied. Daarvoor gebruiken we de library tmap. Daarmee kun je vergelijkbaar met ggplot een kaart samenstellen.

require(tmap)
tm_shape(energielabels) +
  tm_polygons("MEAN_label", style="jenks", alpha=.5, palette=colorRampPalette(c("green", "red"))(5)) +
  tm_compass(type="arrow", position=c("right", "top"), fontsize = 2 ) + 
  tm_scale_bar()

## Kaart Blokverwarming Nu proberen we de file met blokverwarming in te lezen met readOGR.

# blokverwarming <- readOGR(f_blokverwarming)   # geeft een error

Die geeft een error, omdat blokverwarming een multiploint shapefile is. Dat formaat kan niet worden ingelezen met package gdal. We gebruiken st_read uit package sf. (NB readOGR heeft de voorkeur, mits die bruikbaar is natuurlijk.) Vervolgens plotten we de punten met behulp van base plot.

require(sf)
blokverwarming <- st_read(f_blokverwarming)
## Reading layer `blokverwarming_DH' from data source `C:\Users\Elisabeth\Desktop\Fourpoints\Coursera\GIS_in_R_voorbeeld_Den_Haag\blokverwarming_DH.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 2 fields
## geometry type:  MULTIPOINT
## dimension:      XY
## bbox:           xmin: 4.222274 ymin: 52.01797 xmax: 4.398967 ymax: 52.11275
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
plot(blokverwarming$geometry,pch='.')

## Kaart Energielabels met Blokverwarming Met behulp van base plot kunnen we de punten van blokverwarming over de p5 van energielabels heen.

# wat gevogel met x en y coordinaten uit het sf object halen
x_blok <- rep(1.1,length(blokverwarming$geometry[[1]])/2)
y_blok <- rep(1.1,length(blokverwarming$geometry[[1]])/2)
for (i in 1:length(blokverwarming$geometry[[1]])/2)
{
  x_blok[i] <-  blokverwarming$geometry[[1]][i,][1]
  y_blok[i] <-  blokverwarming$geometry[[1]][i,][2]
  i=i+1
}
blokverwarming_xy <- as.data.frame(cbind(x_blok,y_blok))
# plot
plot(energielabels,col=colorRampPalette(c("brown", "grey"))(5))
points(blokverwarming_xy,pch='.',col="blue")

We willen natuurlijk liever zo’n mooie tm_map plot maken.

blokverwarming_sp <- SpatialPoints(coords=blokverwarming_xy, proj4string = CRS(proj4string(energielabels) ) )
tm_shape(energielabels) +
  tm_polygons("MEAN_label", style="jenks", alpha=.5, palette=colorRampPalette(c("green", "red"))(5)) +
  tm_shape(blokverwarming_sp) +
  tm_dots(col="blue",legend.show = FALSE, legend.is.portrait = FALSE ) +
  tm_compass(type="arrow", position=c("right", "top"), fontsize = 2 ) + 
  tm_scale_bar()

## Open Street Map als referentie Alleen de vormen van p5 gebieden zijn eigenlijk niet genoeg voor een kaart. Om een basisreferentie van Den Haag/Nederland te hebben nemen we een OSM als achtergrond. We kunnen bijvoorbeeld ggmap gebruiken.

require(ggmap)
# Den Haag longitude and latitude
DH_loc = c(lon = 4.3, lat = 52.05 )
DH_basemap <- get_map(location = DH_loc, zoom = 12 )
ggmap(DH_basemap)

Maar met de instelling ‘tmap_mode(“view”)’ kun je een in-/uitzoombare kaart maken waar meteen al een basemap ondergeplot wordt. (Alleen kan dan het kompas niet laten zien worden.) Als je helemaal inzoomt kun je precies zien welk kadasterveld blokverwarming heeft.

tmap_mode("view")
tm_shape(energielabels) +
  tm_polygons("MEAN_label", style="jenks", alpha=.8, palette=colorRampPalette(c("green", "red"))(5)) +
  tm_shape(blokverwarming_sp) +
  tm_dots(col="blue",size=0.01) +
  tm_scale_bar()

Aantal panden met blokverwarming per p5

Een veel uitgevoerde bewerking bij het maken van kaarten is spatial join. Hierbij worden twee verschillende kaarten of layers samengevoegd. We voeren nu een spatial join uit van de p5 gebieden met energielabels, en de punten van de blokverwarming. En we maken een kaart met kleurindex op basis van aantal panden met blokverwarming.

a=rep(99999,length(energielabels@polygons))
for ( i in 1:length(energielabels@polygons))
{
  a[i]=sum(point.in.polygon(blokverwarming_xy$x_blok,
                   blokverwarming_xy$y,
                   energielabels@polygons[[i]]@Polygons[[1]]@coords[,1],
                   energielabels@polygons[[i]]@Polygons[[1]]@coords[,2]))
  i=i+1
}
energielabels$N_blokverwarming<-a
tm_shape(energielabels) +
  tm_polygons("N_blokverwarming", style="jenks", alpha=.8, palette=colorRampPalette(c("lightyellow", "brown"))(10)) +
  tm_shape(blokverwarming_sp) +
  tm_dots(col="blue",size=0.01) +
  tm_scale_bar()

Laten we ze naast elkaar plotten om ze te vergelijken.

tmap_mode("plot")
plot1 <- tm_shape(energielabels) +
  tm_polygons("MEAN_label", style="jenks", alpha=.8, palette=colorRampPalette(c("green", "red"))(5))+
  tm_compass(type="arrow", position=c("right", "top"), fontsize = 2 ) +
  tm_scale_bar()
plot2 <- tm_shape(energielabels) +
  tm_polygons("N_blokverwarming", style="jenks", alpha=.8, palette=colorRampPalette(c("lightyellow", "brown"))(10)) +
  tm_compass(type="arrow", position=c("right", "top"), fontsize = 2 ) +
  tm_scale_bar()
require(grid)
grid.newpage()
pushViewport(viewport(layout=grid.layout(1,2)))
print(plot1, vp=viewport(layout.pos.col = 1))
print(plot2, vp=viewport(layout.pos.col = 2))

Het lijkt niet gecorreleerd te zijn, maar we testen het lekker toch. Eerst kijken we of energielabel normaal verdeeld is. Helaas is dat niet zo, ook niet na een logtransformatie.

hist(energielabels$MEAN_label)

hist(log10(energielabels$MEAN_label))

## Raster Satelliet Data Op de website van de VU (“http://geoplaza.vu.nl/data/dataset/topraster”) heb ik een raster kaart gevonden (.tif). Helaas lukt het niet om deze vanuit R te downloaden, maar de file staat in de github. Het gaat om topografische rasterdata. Dit houdt in dat elk vierkantje (raster) een kleur heeft gekregen voor een bepaald soort bodemgebruik. Deze kleur bestaat uit drie waarden (RGB). Het gaat hier dus om een 3 layer .tif. We lezen de drie lagen apart in, en plotten ze. Als je R zelf de kleuren laat kiezen krijg je bij elke laag dezelfde, ‘groenig’.

require(raster)
r1 = raster("Top_250_raster_(2009)-1516618615.tif", band = 1)
r2 = raster("Top_250_raster_(2009)-1516618615.tif", band = 2)
r3 = raster("Top_250_raster_(2009)-1516618615.tif", band = 3)
plot(r1); plot(r2); plot(r3)

Eigenlijk hebben we daar niet zoveel aan. We willen de drie lagen natuurlijk over elkaar bekijken. Met het stack() commando kunnen we drie lagen in één keer inlezen. Met plotRGB() kunnen we de drie lagen over elkaar heen plotten, en krijgen we de keuren te zien zoals ze bedoeld zijn.

r=stack("Top_250_raster_(2009)-1516618615.tif")
plotRGB(r)

Voor de grap kunnen we ook een spatialJoin doen, bijvoorbeeld de gemiddelde waardes (RGB) per stadsdeel berekenen. Gewoon puur om te laten zien dat het kan, want een dergelijke join geeft natuurlijk helemaal geen informatie. (Pas op, duurt lang om dit chunk te runnen…)

# met extract kunnen we per polygon een functie toepassen. bijvoorbeeld mean.
v1 <- extract(r1, energielabels, fun = mean)
v2 <- extract(r2, energielabels, fun = mean)
v3 <- extract(r3, energielabels, fun = mean)
v1 <- round(v1,0)
v2 <- round(v2,0)
v3 <- round(v3,0)
output = data.frame(cbind(v1,v2,v3))
names(output) = c("R","G","B")

R print de waarschuwing: “Transforming SpatialPolygons to the CRS of the Raster”. R heeft namelijk gezien dat de projectie (CRS/proj4string) van beide objecten anders is. Een raster is erg lastig van projectie te veranderen, daarom veranderd R de projectie van de polygons. Het restullaat is zoals verwacht tamelijk nietszeggend. Al lijkt het nog enigszins op het originele raster.

energielabels@data <- cbind(energielabels@data,output)
plot(energielabels,col=rgb(energielabels@data$R,energielabels@data$G,energielabels@data$B,maxColorValue = 255))

## Panden en wegvakken Amsterdamse Veerkade Martijn en René hebben eerder dit jaar twee shapefiles gemaakt van de Amsterdamse Veerkade: panden, of eigenlijk gevels; en wegvakken, of eigenlijk polygonen. Laten we deze eens inlezen en bekijken in R (base plot).

# readOGR op de panden (= lijnen)
pand <- readOGR("pand_data_geo_merged.shp")
## OGR data source with driver: ESRI Shapefile 
## Source: "pand_data_geo_merged.shp", layer: "pand_data_geo_merged"
## with 88 features
## It has 21 fields
## Integer64 fields read as strings:  FID_1 ID Count_
# readOGR op de straten (= polygons)
wegvak <- readOGR("Wegvak_geo_merged.shp")
## OGR data source with driver: ESRI Shapefile 
## Source: "Wegvak_geo_merged.shp", layer: "Wegvak_geo_merged"
## with 44 features
## It has 68 fields
plot(wegvak, col="lightblue")
lines(pand, lty=5, col="red", lwd=3)

En in tmap:

tmap_mode("view")
tm_shape(wegvak) +
  tm_polygons("Avg_PM10_O",style="jenks", palette=colorRampPalette(c("green", "red"))(5)) +
  tm_shape(pand) +
  tm_lines("Avg_AANSL", lwd = 10, lty = "dotted", style="jenks", palette=colorRampPalette(c("green", "red"))(5)) +
  tm_scale_bar()

Er zijn een paar gekke dingetjes aan de hand. Aan de wegvakken met een lage PM10 zitten de panden(gevels) met de laagste aantallen slaapgestoorden. Er kunnen twee dingen aan de hand zijn: óf de data in de shapefile is incorrect, óf de gevels aan de PM10-rijke straten zijn voornamelijk commercieel. ## Work in Progress Geluid

## Er stond ook een .tiff bestand op de dropbox (raster dus)
geluid <- stack("RIS Geluidbelasting Veerkade 2011.tiff")
plotRGB(geluid)

geluid_legenda <- stack("RIS Geluidbelasting Veerkade 2011 Legenda.tiff")
plotRGB(geluid_legenda)